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1. Introduction 


The equations 


V 2 co = 0, 

w = _v 2 x 



(1.1a) 

(1.1b) 


describe, in suitable units, two-dimensional Stokes flow of an incompressible fluid 
occupying a domain D in which to is the vorticity and X is the stream function. The flow is 
uniquely determined by specifying the velocity on the boundary B of D, a condition which 
leads to specifying the stream function X and its normal derivative Xn on B. A 
mathematically similar problem arises in describing the equilibrium of a flat plate in 
structural mechanics where a related one-dimensional problem describes the equilibrium of 
a clamped beam. A key to treating these simple problems by finite difference or finite 
element methods is to introduce effective methods for imposing the boundary conditions 
through which (1. la) is coupled to (1.1b). These models thus provide a simple starting 
point for examining the general treatment of boundary conditions for more general time- 

dependent Navier-Stokes incompressible flows. 
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For the purpose of our discussion we may assume D is a square domain. A 
standard finite difference method to solve ( 1 . 1 ) is to introduce a uniform grid and then 
employ standard five-point finite difference operators to express each equation in (1.1). At 
any point on the boundary B a value of X is specified by the boundary conditions but a 
value of co at the same boundary mesh point will also be required to complete the 
computation. Peyret and Taylor [1] review the use of extrapolation methods to achieve this 
using Taylor series arguments. However, this technique can be expected to be of limited 
value when finte-volume methods are used to treat curved boundaries or when geometrical 
singularities arise when using curvilinear coordinates. As we shall also see, its use can be 
expected to result in a loss of accuracy even in simple cases. The method discussed in this 
report can be expected to overcome these difficultites. 


2. A Compact Difference Scheme 

We first describe a compact finite difference scheme for solving 4>"= g which has 
been described in Rose[2]. It is a specialization to one-dimension of a more general finite 
volume scheme for solving div u = g, grad (J) = u on general domains. 

The scheme expresses a relationship between certain primary variables 4>j, u (xj), 
u + (xj) which are associated with the endpoints of the non-over lapping intervals (x^ Xj+i), 
i = 0,1,...,M-1, and dual variables 4>j which are associated with the interval centerpoints 
Xj, j = 1/2, 3/2,.. .,M- 1/2, which we regard as forming a dual (or staggared) grid. (In the 
following, the index i will be associated with the primary grid, while j will indicate the dual 
grid.) Indicating the cell endpoints by the symbol x and the cell midpoints by the symbol 
o, the points form the pattern 

xoxoxoxoxo X....X o x 
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totalling 2M+1 points. 

For uniform grids the scheme to solve 4> = g can be described as follows: write the 
equation as the first order system u’ = g, 4>' = u and express the first equation by the 
difference equations 

[ u _ (xj+i/2) - u + (xj_i/2)] / Ax = gj j = 1/2, 3/2,..., M-l/2 

(2.1) 

in each interval of length Ax = 2h. We interpret the second equation as relating u~ andu + 
with forward and backward differences involving the variables ({) at primary and secondary 

points of the grid: 

u-(Xj) = (4^ -4 ^_i/ 2 ^/ i = 1,2,..., M (2.2) 

u + (xj) = (4> 1+ 1/2 -<J>,)/h, » = M_1 

Introduce the central average and difference operators 

M 4 $= ( 4 >j+l /2 + 4 ^-l/ 2) /2 > A 4 ^= ( 45+1/2 - 45-1/2)- 

and impose the condition that the variables u~(xj) and u + (xj) be continuous at interior 

endpoints of the primary grid, i.e., 

Uj = u-(xj) = u + (xj), i = 1,2,..., M-l. (2.3) 

Using the definitions (2.2) in (2.1) we find, with K = 1/Ax^, 

gj = 4 k (ju<|>j - 45 ). j = 1/2. 3/2,..., M-l/2 (2.4) 

whike use of the continuity conditions (2.3) leads to 

0 = - 4>j i = 1,2,..., M-l (2.5) 

This tri-diagonal system is to be solved for the variables 4) with specified boundary 
conditions of either Dirichlet or Neumann type. 

A cyclic (odd-even) reduction technique leads to the system of equations 

pg^xAfyj, i = 1 , 2 ,..., M-l (2.6) 
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(2.7) 


for the primary variables and, separately, to 

gj = k , j = 3 / 2 , , 5 / 2 .., M-3/2 

for the dual variables. Under Dirichlet-type boundary conditions, the solution of the first 
set of equations can be solved directly for the primary variables; the values <$>y 2 , 4 > m - 1/2 
for the dual variables can be found in terms of these primary values by solving each of the 
equations 

gj = k , j = 1/2, M-l/2 (2.8) 

and these, in turn, can then be used to provide Dirichlet data to solve (2.7). This reduction 
is less useful when Neumann-type data are imposed; in this case it is best to solve (2.4)- 
(2.5) directly. 

The following table compares numerical results obtaained by solving <f> = g by a 
standard finite difference scheme and by the compact scheme just described. 


Error Norm Comparison of Standard and Compact Schemes 
for the Equation 4>"=g 


<Jj=x(1 -x): 

(The precision of results for this example is questionable because of machine limitations) 


endpoint error norms 


standard 

solution derivative 


# intervals 

6 

12 

24 

48 


2.77556e-17 
5.551 12e-17 
1.1 1022e-16 
2.63678e-16 


.166667 

8.33333e-2 

4.16667e-2 

2.08333e-2 


solution 

5.55 1 12e- 17 
8.32667e-17 
2.49800e-16 
1 .83 1 87e- 1 5 


compact 

derivative 

5.55 1 1 2e- 1 7 
3.33067e-16 
9.99201e-16 
5.7731 6e- 1 5 
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<t> = x 2 (1-x) 2 : 


endooint error norms 

standard compact 


# intervals 

solution 

derivative 

solution 

derivative 

6 

6.17284e-3 

6.48148e-2 

1.23457e-2 

1.85185e-2 

12 

1.7361 le-3 

4.16667e-2 

3.47222e-3 

9.25926e-3 

24 

4.34028e-4 

2.40162e-2 

8.68056e-4 

2.89352e-3 

48 

1.08507e-4 

1.62399e-2 

2.17014e-4 

7.957 18e-4 


3. The Clamped Beam Problem 


The deflection <Kx) of a uniform, straight beam under a load -fl(x) per unit length on an 
interval [1_, 1+] is, in dimensionless form, governed by the simple fourth order differential 
equation 

(JT" = -f (3.1) 

If the slope u, bending moment v, and shear force w are given by 

u = 4>' , v = -u' , w = v' (3.2) 

then typical well-posed boundary conditions allow one to prescribe pairs of values among 
(4>,u,v,w) at each endpoint of the interval. We may also write (3.1) as a coupled system of 
second-order equations for (J> and v 

(a) v" = f, (3-3) 

(b) <|>" = -v. 

For a problem in which the boundary conditions involve the pair of values (4), v) at each 
endpoint, for example, a simple Green’s function technique allows (a) to be solved for v(x) 
in terms of f(x) together with the boundary conditions v(l±), while a similar construction 
gives <Kx) in terms of the boundary conditions 4>(1±) and v(x). However, for clamped 
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boundary conditions, 4> = u = 0 at both endpoints and this Green’s function construction 
fails. 


A standard finite difference technique for handling clamped boundary conditions 
involves the use of a Taylor series expansion at each endpoint in order to express v(l±) in 
terms of 4/l±) and u(l±), both values of which are prescribed, as well as one or more 
values of 4> at interior mesh points of the interval. We shall describe this technique first, 
using a rather direct finite difference argument which is suggested by a method described 
in Peyret and Taylor [ 1 ] in connection with a treatment of the Navier-Stokes equations in 
vorticity-stream function variables. We then discuss and illustrate another, closely related, 
technique which arises from an application of a compact finite difference method (Rose[2]) 
to this problem. 

4. Some Standard Finite Difference Approaches . 

We adopt the standard finite-difference notations Xj = i Ax, h = Ax/2, u(xj) = Uj. 
Divide [1_, 1 + ] into M non-overlapping intervals Ij = {x| Xj_i /2 ^ x ^ Xj+i^} with 
centerpoints xj, j= 1/2, 3/2,..., M-l/2. Also, recall the central average and difference 
operators 

P 4 >j= ( <t>j+i /2 + 4 >j-l/ 2 V 2 > A <t>j= ( 45+1/2 - 45-1/2)- ( 4 -!) 

introduced earlier. 

With K = 1/Ax 2 ,a standard finite difference approach to solving (4.3) is to consider 
the coupled difference equations 


(a) fj = k A 2 Vj i = 1, 2, ..., M-l (4.2) 

(b) -Vj = k 
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each of which may be separately solved with Dirichlet-type data. For the clamped beam, 
the boundary conditions <h(l±) = 0, u(l±) = 0 translate into 

(a) 4\) = 4 >m = 0, (4.3) 

(b) uq = u M = 0. 

Were Vj known, (4.2b) could easily be solved for 4) under either pair of these boundary 
conditions by a standard tri-diagonal solver (in case (b) we can add the condition <ho = 0) 

In order to solve (4.2a) we will specify values of v 0 and v M Consistency requires that 
these values be related to values of 4> and u at points on or near the boundary points and we 
may write 

v 0 = B 0 (4>,u), (4.4) 

V M = B m (({) ’ u) 

where B is a suitable boundary operator which incorporates the prescribed boundary 
conditions for.<t> and u. 

One method to obtain a boundary operator B is to use a Taylor series approximation 
as follows: Write 

<J >1 =4> 0 + Ax u o - MAx) 2 Vo +... (4.5) 

4>M-l =( t>M _Axu M- MAx) 2 v m +..., 

so that by imposing the homogeneous clamped beam boundary conditions on 4> and u we 
find 

2 (hj = - (Ax) 2 vq + ... (4-6) 

2< t>M-l = " (Ax) 2 vm +.... 

Higher order extrapolations are discussed in Peyret and Taylor[ 1 ]. 
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5. A Compact Scheme for the Clamped Beam . 


Our objective will be to describe a difference scheme based upon these ideas which 
solves (1) under the clamped beam boundary conditions stated in (4.3), i.e., 

4)0 = 4^ = 0, uq = u M = 0. The boundary conditions for v will not require the additional 
use of a Tavlor series other than that which is implicit in the proposed differe nce equations. 

We will consider the coupled system 

-jLtVj = k , i = 1,2,..., M-l (5.1) 

fj = k A 2 Vj , j = 3/2, ,5/2.., M-3/2. 

The first of these corresponds to (4.5) for the primary variables associated with a compact 
scheme for solving 4> " = - v and the second corresponds to (4.6) for the dual variables 
associated with a compact scheme for solving v" = f. 

Using the definitions of u given by (4.2) the values of v at the endpoints of the dual 
grid are found to satisfy 

- Ax Vj = (uj+i /2 - vij_ 1 / 2 ) j = 1/2 » M-l/2 (5.2) 

and imposing the Neumann-type boundary conditions for u gives 

uj = - Ax v 1/2 , u M _! = Ax v M _ 1/2 . 

Using the definition of u leads to 

2k (4>j -4>i/2)= -V 1/2 (5.3) 

2K ^M-l/2 “ V M-l/2' 

One difference between this and the treatment of boundaiy conditions for the 
standard finite difference method described earlier lies in the position on the mesh of at 
which the boundary values assigned to v are imposed. In the standard problem Vj was 
considered a primaiy variable and a Taylor expansion was required in order to furnish an 
additional equation which allowed boundary values of v at the endpoints of the domain to 
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be determined by cfy. In the compact scheme, on the other hand, the Vj are dual variables 
and the boundary conditions , which determine their values at the boundary of the dual 
mesh, are a consequence of applying the difference equations jlivj = - k A 2 ^ when 
i = 1/2, M-l/2. The cyclic reduction technique allowed the dual variables to be solved 
directly with such data. 

This distinction between variables defined on a primary and dual grid is a common 
feature in treating the Navier-Stokes equations in primitive variables where the pressure 
term is commonly associated with the dual grid. This connection is worth further 
exploration. 

Numerical Examples for the Clamped Beam . 

The following table compares numerical results for the clamped beam problem 
using the standard finite difference method with extrapolation techniques to set the 
boundary conditions and the compact scheme described above. 

Error Norm Comparison of Standard and Compact Schemes 
for a Clamped Beam 

test solution : (J)=x2(l-x)2: 

Standard scheme using 1st order extrapolated BC; 


# intervals 

f 

u 

w 

12 

1.07692 

.213675 

1.92901e-2 

24 

.324713 

7.24377e-2 

1.17071e-2 

48 

.124513 

5.87587e-2 

7.65804e-3 

96 

5.41 195e-2 

3.55117e-2 

4.37644e-3 
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Standard scheme using 2nd order extrapolated BC; 


# intervals 

f 

u 

w 

12 

2.80765e-2 

.14098 

.863799 

24 

1.77672e-2 

9.04946e-2 

.313276 

48 

9.78894e-3 

6.49779e-2 

.1223 

96 

5.07409e-3 

3.71921e-2 

5.17336e-2 


ComDact scheme : 




# intervals 

f 

u 

w 

12 

2.82377e-2 

3.8407e-5 

6.9471 le-3 

24 

6.9977e-3 

2.25417e-6 

1.73624e-3 

48 

1.74305e-3 

1.58858e-7 

4.34032e-4 

96 

4.3588 le-4* 

1.399 l2e-7* 

1.08463e-4 


♦precision doubtful because of machine limitations. 


Note that, as predicted, the compact scheme furnishes second-order accuracy for 
the variables. 


6. A Time-dependent Stokes-tvoe Problem 


The time-dependent equations 

w t = V 2 to, 


to = -v 2 x 


can serve as a model for studying the effect of handling boundary conditions and differs 
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from the Navier-Stokes equations in that only the convective terms have been omitted. A 
method for adapting ADI techniques to solve this problem by a compact scheme is outlined 
in [2]. 

A code to test the accuracy of solving this problem by a compact scheme has been 
developed by J. M. Klimkowski and is reported upon separately. 
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